//==============================================================================
// SequenceAlignment.cc
//
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "SequenceAlignment.hh"
#include <iomanip.h>

static long limitingCaseCount=0;

// KJD Note: This should be changed to 0,1000000, and 2000000
// after testing.
#define STable 0
#define XTable 10
#define YTable 20
#define NINFINITY -10000

void StartTiming(int &a);
int  EndTiming(int b);


//---------------------------------------------------
// * StartTiming
//---------------------------------------------------
//
void 
SequenceAlignment::StartTiming(int &StartTime)
{
  StartTime = clock();
}

//---------------------------------------------------
// * EndTiming
//---------------------------------------------------
//
int 
SequenceAlignment::EndTiming(int StartTime)
{
  long endTime = clock();
  long elapsedTime = (long)(endTime - (long)StartTime);
  int seconds =  (int)((float)elapsedTime/(float)CLOCKS_PER_SEC);
  return(seconds);
}

//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
void
SequenceAlignment::SetScoreParams(int match,
						 int mismatch,
						 int gap)
{
	mMatchScore = match;
	mMismatchScore = mismatch;
	mGapScore = gap;
}

//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
void
SequenceAlignment::SetScoreParams(int match, int mismatch, int gap, int extension)
{
	mMatchScore = match;
	mMismatchScore = mismatch;
	mGapScore = gap;
	mGapExtensionScore = extension;
}


//set
void SequenceAlignment::SetBandScore(int score)
{
	mBandScore=score;
}

//if it is set to 1, then N is treated as N, otherwise, N is treated
//as any letter, always match. Default is set at 1
void SequenceAlignment::SetIgnorN(int choice)
{
	mIgnorN=choice;
}

int SequenceAlignment::IgnorN()
{
	return mIgnorN;
}

//----------------------------------------
//SequenceAlignment constructor
//---------------------------------------
//
SequenceAlignment::SequenceAlignment()
{
  SetScoreParams(MATCH_SCORE,MISMATCH_SCORE,GAP_SCORE);
  SetIgnorN(1);
  SetBandScore(-100000);
}


//------------------------------------------------------
// * SequenceAlignment
//------------------------------------------------------
// 
SequenceAlignment::SequenceAlignment(Sequence S,Sequence T)
{
  mS.clear();
  mT.clear();
  mS = S;
  mT = T;
  mRows = S.size()+1;
  mCols = T.size()+1;
  SetScoreParams(MATCH_SCORE,MISMATCH_SCORE,GAP_SCORE);
  SetIgnorN(1);
  SetBandScore(-100000);
}

//------------------------------------------------------
// * SequenceAlignment
//------------------------------------------------------
// 
SequenceAlignment::SequenceAlignment(string S,string T)
{
	mS.InitFromString(S);
	mT.InitFromString(T);
	mRows = mS.size()+1;
	mCols = mT.size()+1;
	SetScoreParams(MATCH_SCORE,MISMATCH_SCORE,GAP_SCORE);
	SetIgnorN(1);
	SetBandScore(-100000);

	mScoreTable.Allocate(mRows,mCols);
	mPtrTable.Allocate(mRows,mCols);
}

//------------------------------------------------------
// * ~SequenceAlignment
//------------------------------------------------------
// 
SequenceAlignment::~SequenceAlignment(){
}

//------------------------------------------------------
// * GetGlobalAlignment
//------------------------------------------------------
// Globally aligns two sequences.  
//
// This version has constant gap penalties, so it doesn't distinguish 
// between:
//
//  caggggactacctggaaaataatac
//  ----ggact---tgg------a-a-
//
// and
//
//  caggggactacctggaaaataatac
//  ----ggact---tggaa--------
//
// An affine gap penalty (among others) would prefer the latter. 
//
void 
SequenceAlignment::GetGlobalAlignment(Sequence &AlignS,Sequence &AlignT, int &score)
{
  ComputeGlobalAlignmentScore();
  RecoverGlobalAlignment();
  AlignS = mAlignS;
  AlignT = mAlignT;
  score = mScoreTable[mRows-1][mCols-1];
}

//------------------------------------------------------
// * FindWithLocalAlignment
//------------------------------------------------------
// Use local alignment to determine if the target sequence
// is contained in the source sequence.  This is the method to be called
// Source is the database to search, which is mT,  and Target is the query which give to mS
int
SequenceAlignment::GetLocalAlignment(Sequence &Source, Sequence &Target, int &score)
{	
  // Ensure that they are initially clear. 
  mS.clear();
  mT.clear();
  mS = Source;
  mT = Target;
  mRows = Source.size()+1;  //one row and col longer than the actual sequence length
  mCols = Target.size()+1;
  mAlignStart = -1;
  mMaxScore = 0;

  // Ensure case is same.
  mS.ToLowerCase();
  mT.ToLowerCase();
  
  int StartTime,EndTime;
  
  // This is done here since FindWithLocalAlignment is called with the 
  // source and target directly. 
  mScoreTable.Allocate(mRows,mCols);
  mPtrTable.Allocate(mRows,mCols);
  
  ComputeLocalAlignmentScore();
  RecoverLocalAlignment();
	
  score = mMaxScore;
  if (score > Source.size()){
    cout<<"Score error. score is "<< score <<endl;
  }

  //mScoreTable.DeAllocate();
  //mPtrTable.DeAllocate();

  return(mAlignStart);
}

//-----------------------------------------------------------
//GetInternalAlignment
// find case that like ---------------------------------------
//                              --------------
// the stick out of the second sequence will be penalized and count as a gap
//--------------------------------------------------------------------------
void SequenceAlignment::GetInternalAlignment(Sequence &Source, Sequence &Target)
{	
  // Ensure that they are initially clear. 
  mS.clear();
  mT.clear();
  mS = Source;
  mT = Target;
  mRows = Source.size()+1;  //one row and col longer than the actual sequence length
  mCols = Target.size()+1;
  mAlignStart = -1;
  mMaxScore = 0;

  // Ensure case is same.
  mS.ToLowerCase();
  mT.ToLowerCase();
  
  int StartTime,EndTime;
  
  // This is done here since Alignment is called with the 
  // source and target directly.
  mScoreTable.DeAllocate();
  mPtrTable.DeAllocate(); 
  mScoreTable.Allocate(mRows,mCols);
  mPtrTable.Allocate(mRows,mCols);
  
  ComputeInternalAlignmentScore();
  RecoverInternalAlignment();
       

  return;
}


//-----------------------------------------------------------
//GetEndAlignment
// find case that like ---------------
//                              --------------
// It will start from the internal of the first and end at the sencod seq
//--------------------------------------------------------------------------
int SequenceAlignment::GetEndAlignment(Sequence &Source, Sequence &Target)
{	
  // Ensure that they are initially clear. 
  mS.clear();
  mT.clear();
  mS = Source;
  mT = Target;
  mRows = Source.size()+1;  //one row and col longer than the actual sequence length
  mCols = Target.size()+1;
  mAlignStart = -1;
  mMaxScore = 0;

  // Ensure case is same.
  mS.ToLowerCase();
  mT.ToLowerCase();
  
  int StartTime,EndTime;
  
  // This is done here since Alignment is called with the 
  // source and target directly.
  mScoreTable.DeAllocate();
  mPtrTable.DeAllocate(); 
  mScoreTable.Allocate(mRows,mCols);
  mPtrTable.Allocate(mRows,mCols);
  
  ComputeEndAlignmentScore();
  if( RecoverEndAlignment() ) {
    return 1; //if overlap at the end
  }

  return 0;
}

//-----------------------------------------
//
//---------------------------------------------
void SequenceAlignment::GetBandInternalAlignment(Sequence &Source, Sequence &Target, int targetstart, int bandsize)
{	
  // Ensure that they are initially clear. 
  mS.clear();
  mT.clear();
  mS = Source;
  mT = Target;
  mRows = Source.size()+1;  //one row and col longer than the actual sequence length
  mCols = Target.size()+1;
  mbandwidth=bandsize;
  //decide the relationship between them and allocate score table and ptr table
  mTBandStart=targetstart-bandsize;
  mTBandEnd=targetstart+bandsize; 
  if( (mTBandEnd+mS.size()-1)<=mT.size() ) { //the end is inside
      mBRows=mS.size()+1;
      mBCols=mbandwidth*2+3;
  } else if( (mT.size()-mTBandStart+1)<=mS.size()) { //complete query stick out
    mBRows=mT.size()-mTBandStart+2;
    mBCols=mbandwidth*2+3;
  } else {
    mBRows=mS.size() + 1;
    mBCols=mbandwidth*2+3;
  }

  mAlignStart = -1;
  mMaxScore = 0;
  // Ensure case is same.
  mS.ToLowerCase();
  mT.ToLowerCase();
  
  int StartTime,EndTime;
  
  // This is done here since Alignment is called with the 
  // source and target directly.
  mbandScoreTable.DeAllocate();
  mbandPtrTable.DeAllocate(); 
  mbandScoreTable.Allocate(mBRows,mBCols);
  mbandPtrTable.Allocate(mBRows,mBCols);
  
  if(ComputeBandInternalAlignmentScore()) {
    RecoverBandInternalAlignment();
  }
}

int SequenceAlignment::GetBandEndAlignment(Sequence &Source, Sequence &Target, int targetstart, int bandsize)
{	
  // Ensure that they are initially clear. 
  mS.clear();
  mT.clear();
  mS = Source;
  mT = Target;
  mRows = Source.size()+1;  //one row and col longer than the actual sequence length
  mCols = Target.size()+1;
  mbandwidth=bandsize;
  //decide the relationship between them and allocate score table and ptr table
  mTBandStart=targetstart-bandsize;
  mTBandEnd=targetstart+bandsize; 
  if( (mTBandEnd+mS.size()-1)<=mT.size() ) { //the end is inside, so skip if both side internal
    if(targetstart<0) {
      mBRows=mS.size()+1;
      mBCols=mbandwidth*2+3;
    } else {
      cout<<"the query is complete inside of the target"<<endl;
      return 0;
    }
  } else if( (mT.size()-mTBandStart+1)<=mS.size()) { //complete query stick out
    mBRows=mT.size()-mTBandStart+2;
    mBCols=mbandwidth*2+3;
  } else {
    mBRows=mS.size() + 1;
    mBCols=mbandwidth*2+3;
  }

  mAlignStart = -1;
  mMaxScore = 0;
  // Ensure case is same.
  mS.ToLowerCase();
  mT.ToLowerCase();
  
  int StartTime,EndTime;
  
  // This is done here since Alignment is called with the 
  // source and target directly.
  mbandScoreTable.DeAllocate();
  mbandPtrTable.DeAllocate(); 
  mbandScoreTable.Allocate(mBRows,mBCols);
  mbandPtrTable.Allocate(mBRows,mBCols);
  
  if(ComputeBandEndAlignmentScore() && RecoverBandEndAlignment()) {
    return 1;
  }
  return 0;
}

//GetBandOverlapAlignment
//this is either end or internal alignment
//to check if the alignment is good, rely on maxscore
int SequenceAlignment::GetBandOverlapAlignment(Sequence &Source, Sequence &Target, int targetstart, int bandsize)
{	
  // Ensure that they are initially clear. 
  mS.clear();
  mT.clear();
  mS = Source;
  mT = Target;
  mRows = Source.size()+1;  //one row and col longer than the actual sequence length
  mCols = Target.size()+1;
  mbandwidth=bandsize;
  //decide the relationship between them and allocate score table and ptr table
  mTBandStart=targetstart-bandsize;
  mTBandEnd=targetstart+bandsize; 
  if( (mTBandEnd+mS.size()-1)<=mT.size() ) { //the end is inside, so skip if both side internal
      mBRows=mS.size()+1;
      mBCols=mbandwidth*2+3;
  } else if( (mT.size()-mTBandStart+1)<=mS.size()) { //complete query stick out
    mBRows=mT.size()-mTBandStart+2;
    mBCols=mbandwidth*2+3;
  } else {
    mBRows=mS.size() + 1;
    mBCols=mbandwidth*2+3;
  }

  mAlignStart = -1;
  mMaxScore = 0;
  // Ensure case is same.
  mS.ToLowerCase();
  mT.ToLowerCase();
  
  int StartTime,EndTime;
  
  // This is done here since Alignment is called with the 
  // source and target directly.
  mbandScoreTable.DeAllocate();
  mbandPtrTable.DeAllocate(); 
  mbandScoreTable.Allocate(mBRows,mBCols);
  mbandPtrTable.Allocate(mBRows,mBCols);
  
  if(ComputeBandEndAlignmentScore()) {
    RecoverBandEndAlignment();
    return 1;
  } else {
    //cout<<"stop the band"<<endl;
    return 0;
  }
}

#if 0
//------------------------------------------------------
// * Score
//------------------------------------------------------ 
//
inline int 
SequenceAlignment::Score(char A,char B)
{
	int score;
	if (A == B){
		score=mMatchScore;
	}else{
		score=mMismatchScore;
	}
	return(score);
}
#endif


int SequenceAlignment::ComputeBandInternalAlignmentScore() {

	int diagScore,leftScore,topScore, temp;
	int i,j;
	int vrow, vcol;
	char *pS,*pT;
  
	pS = new char[mRows];
	pT = new char[mCols];
	for(i = 0;i<mRows;i++){
	  pS[i]=mS[i];
	}
	for(j=0;j<mCols;j++){
	  pT[j]=mT[j];
	}

	//there is negative infinity penalty outside the band
	for(i = 0;i<mBRows;i++) {
	  mbandScoreTable[i][0] = NINFINITY;
	  mbandScoreTable[i][mBCols-1]= NINFINITY;
	}
	
	for(j=1; j<mBCols-1; ++j) {
	  mbandScoreTable[0][j]=0; //can start anywhere in the target
	}

	//initialize the table so that for all virtual cols less or equal to 0 will be gap penalized
	//for virtal cols greater than 0 and vrow is 0, just with 0 so that match can start anywhere in the target
	for(i=0; i<mBRows; ++i) {
	  vcol=i+mTBandStart-2; //the NIFINITY POS, start of a row
	  if(vcol>0) break;  //end of initilization
	  for(j=1; j<mBCols-1; ++j) {
	    if( (j+vcol)<1) { //initilize upto real vcol is 0
	      mbandScoreTable[i][j]=i*mGapScore;
	    }
	  }
	}

	//the calculation begin, instead of left, top and diagnal, it is left, righttop and top
	for(i=1;i<mBRows;i++){
	  vcol=i+mTBandStart-2;
	  for(j=1;j<mBCols-1;j++){
	    if( (j+vcol)<1 ) continue; //skip those initilized ones
	    if( (j+vcol) > mT.size() ) continue; //skip ones that outside target boundary
	    diagScore = mbandScoreTable[i-1][j]+ Score(pS[i-1],pT[vcol+j-1]);
	    leftScore = mbandScoreTable[i][j-1]  + mGapScore;
	    topScore =  mbandScoreTable[i-1][j+1]  + mGapScore;
			
	    temp=Max(diagScore,leftScore, topScore);
	    if(temp<mBandScore) return 0; //if the score is way to low, stop right away
	    mbandScoreTable[i][j]=temp;
	    // Save pointers to the places where this value could have
	    // come from.
	    mbandPtrTable[i][j]=0; //clear the ptrtable
	    if (temp == diagScore) {
	      Set(DIAG,mbandPtrTable[i][j]);
	    } else if (temp == leftScore) {
	      Set(LEFT,mbandPtrTable[i][j]);
	    } else if (temp == topScore) {
	      Set(TOP,mbandPtrTable[i][j]);
	    }
	  }
	}
	delete [] pS;
	delete [] pT;
	return 1;
	//PrintBandScoreTable();
}

//ComputeBandEndAlignmentScore
//

int SequenceAlignment::ComputeBandEndAlignmentScore() {

	int diagScore,leftScore,topScore, temp;
	int i,j;
	int vrow, vcol;
	char *pS,*pT;
  
	pS = new char[mRows];
	pT = new char[mCols];
	for(i = 0;i<mRows;i++){
	  pS[i]=mS[i];
	}
	for(j=0;j<mCols;j++){
	  pT[j]=mT[j];
	}

	//there is negative infinity penalty outside the band
	for(i = 0;i<mBRows;i++) {
	  mbandScoreTable[i][0] = NINFINITY;
	  mbandScoreTable[i][mBCols-1]= NINFINITY;
	}
	
	for(j=1; j<mBCols-1; ++j) {
	  mbandScoreTable[0][j]=0; //can start anywhere in the target
	}

	//initialize the table so that for all virtual cols less or equal to 0 will be 0
	//for virtal cols greater than 0 and vrow is 0, just with 0 so that match can start 
	//anywhere in the target and query
	for(i=0; i<mBRows; ++i) {
	  vcol=i+mTBandStart-2; //the NIFINITY POS, start of a row
	  if(vcol>0) break;  //end of initilization
	  for(j=1; j<mBCols-1; ++j) {
	    if( (j+vcol)<1) { //initilize upto real vcol is 0
	      mbandScoreTable[i][j]=0;
	    }
	  }
	}

	//the calculation begin, instead of left, top and diagnal, it is left, righttop and top
	for(i=1;i<mBRows;i++){
	  vcol=i+mTBandStart-2;
	  for(j=1;j<mBCols-1;j++){
	    if( (j+vcol)<1 ) continue; //skip those initilized ones
	    if( (j+vcol) > mT.size() ) continue; //skip ones that outside target boundary
	    diagScore = mbandScoreTable[i-1][j]+ Score(pS[i-1],pT[vcol+j-1]);
	    leftScore = mbandScoreTable[i][j-1]  + mGapScore;
	    topScore =  mbandScoreTable[i-1][j+1]  + mGapScore;
			
	    temp=Max(diagScore,leftScore, topScore);
	    if(temp<mBandScore) return 0;
	    mbandScoreTable[i][j]=temp;
	    // Save pointers to the places where this value could have
	    // come from.
	    mbandPtrTable[i][j]=0; //clear the ptrtable
	    if (temp == diagScore) {
	      Set(DIAG,mbandPtrTable[i][j]);
	    } else if (temp == leftScore) {
	      Set(LEFT,mbandPtrTable[i][j]);
	    } else if (temp == topScore) {
	      Set(TOP,mbandPtrTable[i][j]);
	    }
	  }
	}
	delete [] pS;
	delete [] pT;
	//PrintBandScoreTable();
	return 1;
}

//------------------------------------------------------
// * ComputeLocalAlignmentScore
//------------------------------------------------------
// Optimized version of ComputeAlignmentScore.  Clarity is
// sacrificed, a little bit, for the sake of speed.  
//
//
void 
SequenceAlignment::ComputeLocalAlignmentScore()
{
	register int MaxScore = 0,MaxRow=0,MaxCol=0;
	register int i,j,im,jm,rows,cols,temp,diagScore,leftScore,topScore;
	register char *pS,*pT;
	
	// Native character array access is MUCH faster than vector<char>
	// access.  To illustrate, this single efficiency measure dropped
	// the running time from 8 sec to 3 sec, suggesting that 60% of 
	// the inner loop time was spent accessing these characters. 
	pS = new char[mRows];
	pT = new char[mCols];
	for(i = 0;i<mRows;i++){
	  pS[i]=mS[i];
	}

	for(j=0;j<mCols;j++){
	  pT[j]=mT[j];
	}
	
	rows = mRows;
	cols = mCols;

	// Init first two rows to zero for local alignment
	for(i = 0;i<rows;i++)
		mScoreTable[i][0] = 0;
	for(j = 0;j<cols;j++)
		mScoreTable[0][j]=0;

	//-------------- Speed critical section ---------------------------
	//
	// Everything in this loop will happen tens of millions of times
	// so it is critical that the total number of operations be minimized.
	// This justifies the expense of copying everything to local register
	// variables, etc.  I believe that this can be optimized further by 
	// combining Max and the subsequent Set functions, and by ordering 
	// the comparisons to take advantage of the fact that the score
	// in a local alignment is most everywhere zero.
	//
	for(i=1;i<rows;i++){
	  for(j=1;j<cols;j++){
		  
		 im = i-1;
		 jm = j-1;
		 diagScore = mScoreTable[im][jm] + Score(pS[im],pT[jm]);
		 leftScore = mScoreTable[i][jm]  + mGapScore;
		 topScore =  mScoreTable[im][j]  + mGapScore;
		 
		 temp = Max(topScore,leftScore,diagScore,0);

		 // KJD BUG HERE??? Note that temp can be zero at times when 
		 // temp != diagScore, leftScore, OR topScore.  What then???
		 // The pointer in the table gets set to nothing in particular. 
		 // What should happen??
	 
		 mScoreTable[i][j]=temp;

		 mPtrTable[i][j]=0;
		 if (temp == topScore) Set(TOP,mPtrTable[i][j]);
		 if (temp == leftScore) Set(LEFT,mPtrTable[i][j]);
		 if (temp == diagScore) Set(DIAG,mPtrTable[i][j]);
		 
		 // Keep track of the location of the maximum score. 
		 if (temp > MaxScore){
			MaxScore = temp;
			MaxRow = i;
			MaxCol = j;
		 }
	  }
	}
	//-------------- End speed critical section -----------------------
	
	mMaxScore = MaxScore;
	mMaxRow = MaxRow;
	mMaxCol = MaxCol;
	delete [] pS;
	delete [] pT;
}

//------------------------------------------------------
// * ComputeIntrenalAlignmentScore
//------------------------------------------------------
// 
void 
SequenceAlignment::ComputeInternalAlignmentScore()
{
	int diagScore,leftScore,topScore, temp;
	int i,j;
	char *pS,*pT;
	
	// Native character array access is MUCH faster than vector<char>
	// access.  To illustrate, this single efficiency measure dropped
	// the running time from 8 sec to 3 sec, suggesting that 60% of 
	// the inner loop time was spent accessing these characters. 
	pS = new char[mRows];
	pT = new char[mCols];
	for(i = 0;i<mRows;i++){
	  pS[i]=mS[i];
	}

	for(j=0;j<mCols;j++){
	  pT[j]=mT[j];
	}

	//there is gap penalty for start in the middle of the query seq
	for(i = 0;i<mRows;i++)
		mScoreTable[i][0] = i*mGapScore;
	//the target sequence can start from anywhere
	for(j = 0;j<mCols;j++)
		mScoreTable[0][j]= 0;

	for(i=1;i<mRows;i++){
		for(j=1;j<mCols;j++){
			
			diagScore = mScoreTable[i-1][j-1]+ Score(pS[i-1],pT[j-1]);
			leftScore = mScoreTable[i][j-1]  + mGapScore;
			topScore =  mScoreTable[i-1][j]  + mGapScore;
			
			temp=Max(diagScore,leftScore, topScore);
			mScoreTable[i][j]=temp;

			// Save pointers to the places where this value could have
			// come from.
			mPtrTable[i][j]=0; //clear the ptrtable
			if (temp == diagScore) {
			  Set(DIAG,mPtrTable[i][j]);
			} else if (temp == leftScore) {
			  Set(LEFT,mPtrTable[i][j]);
			} else if (temp == topScore) {
			  Set(TOP,mPtrTable[i][j]);
			}
		}
	}
	delete [] pS;
	delete [] pT;
}
	
//------------------------------------------------------
// * ComputeEndAlignmentScore, an exact overlap algorithm
//------------------------------------------------------
// 
void 
SequenceAlignment::ComputeEndAlignmentScore()
{
	int diagScore,leftScore,topScore, temp;
	int i,j;
	char *pS,*pT;
	
	// Native character array access is MUCH faster than vector<char>
	// access.  To illustrate, this single efficiency measure dropped
	// the running time from 8 sec to 3 sec, suggesting that 60% of 
	// the inner loop time was spent accessing these characters. 
	pS = new char[mRows];
	pT = new char[mCols];
	for(i = 0;i<mRows;i++){
	  pS[i]=mS[i];
	}

	for(j=0;j<mCols;j++){
	  pT[j]=mT[j];
	}

	//there is gap penalty for start in the middle of the query seq
	for(i = 0;i<mRows;i++)
		mScoreTable[i][0] = 0;
	//the target sequence can start from anywhere
	for(j = 0;j<mCols;j++)
		mScoreTable[0][j]= 0;

	for(i=1;i<mRows;i++){
		for(j=1;j<mCols;j++){
			
			diagScore = mScoreTable[i-1][j-1]+ Score(pS[i-1],pT[j-1]);
			leftScore = mScoreTable[i][j-1]  + mGapScore;
			topScore =  mScoreTable[i-1][j]  + mGapScore;
			
			temp=Max(diagScore,leftScore, topScore);
			mScoreTable[i][j]=temp;

			// Save pointers to the places where this value could have
			// come from.
			mPtrTable[i][j]=0; //clear the ptrtable
			if (temp == diagScore) {
			  Set(DIAG,mPtrTable[i][j]);
			} else if (temp == leftScore) {
			  Set(LEFT,mPtrTable[i][j]);
			} else if (temp == topScore) {
			  Set(TOP,mPtrTable[i][j]);
			}
		}
	}
	delete [] pS;
	delete [] pT;
}
	
//------------------------------------------------------
// * ComputeGlobalAlignmentScore
//------------------------------------------------------
// 
void 
SequenceAlignment::ComputeGlobalAlignmentScore()
{
	int diagScore,leftScore,topScore,temp;

	mScoreTable.Allocate(mRows,mCols);
	mPtrTable.Allocate(mRows,mCols);
	int i,j;
	for(i = 0;i<mRows;i++)
		mScoreTable[i][0] = i*mGapScore;
	for(j = 0;j<mCols;j++)
		mScoreTable[0][j]=j*mGapScore;

	for(i=1;i<mRows;i++){
		for(j=1;j<mCols;j++){
			
			diagScore = mScoreTable[i-1][j-1]+ Score((char)mS[i-1],(char) mT[j-1]);
			leftScore = mScoreTable[i][j-1]  + mGapScore;
			topScore =  mScoreTable[i-1][j]  + mGapScore;
			
			temp=Max(diagScore,leftScore, topScore);
			mScoreTable[i][j]= temp;
			
			// Save pointers to the places where this value could have
			// come from. 
			if (temp == diagScore) Set(DIAG,mPtrTable[i][j]);
			if (temp == leftScore) Set(LEFT,mPtrTable[i][j]);
			if (temp == topScore) Set(TOP,mPtrTable[i][j]);
		}
	}
}


//------------------------------------------------------
// * Swap
//------------------------------------------------------
void
SequenceAlignment::Swap(int &a,int &b)
{
  int temp = a;
  a = b;
  b = temp;
}


//------------------------------------------------------
// * RecoverGlobalAlignment
//------------------------------------------------------
// Reconstruction of alignment from pointers. 
//  
void 
SequenceAlignment::RecoverGlobalAlignment()
{
  int r = mRows-1,c =mCols-1;

  //reset gap count
  mSGapCount=0;
  mTGapCount=0;
  mMismatchCount=0;
  mAlignS.clear();
  mAlignT.clear();

  // Reconstruct the reverse, then flip it. 
  while((r >0) && (c >0)){
	 if (IsSet(DIAG,mPtrTable[r][c])){
		mAlignS.push_back(mS[r-1]);
		mAlignT.push_back(mT[c-1]);
		if(mS[r-1] != mT[c-1]) mMismatchCount++;
		c--;
		r--;
	 }else{
		if (IsSet(TOP,mPtrTable[r][c])){
		  mAlignS.push_back(mS[r-1]);
		  mAlignT.push_back(GAP);
		  r--;
		  mTGapCount++;
		}else{
		  if (IsSet(LEFT,mPtrTable[r][c])){
			 mAlignS.push_back(GAP);
			 mAlignT.push_back(mT[c-1]);
			 c--;
			 mSGapCount++;
		  }
		}
	 }
  }
	
  // Since we fill back to front, when strings are 
  // substantially different size the beginning will be
  // filled with blanks.  This fills in the alignment for 
  // this section. (only one of r or c should be > 0). 
  if (r > 0){
	 for(int i = r;i> 0;i--) {
		mAlignS.push_back(mS[i-1]);
		mAlignT.push_back(GAP);
	 }
  }
	
  if (c > 0){
	 for(int i = c;i> 0;i--) {
		mAlignS.push_back(GAP);
		mAlignT.push_back(mT[i-1]);
	 }
  }

  // Reverse, since we constructed a backward version.  
  mAlignS.Reverse();
  mAlignT.Reverse();

  //update the start and end of the alignment which are the start and
  //end of the query and target
  mAlignStart=1;
  mAlignEnd=mS.size();
  mTAlignStart=1;
  mTAlignEnd=mT.size();
}

//------------------------------------------------------
// * RecoverInternalAlignment
//------------------------------------------------------
// Reconstruction of alignment from pointers. 
// only accept the alignment that end at the end of the second
// seq, that is, the row=mRow-1
// 
void 
SequenceAlignment::RecoverInternalAlignment()
{
  int r = mRows-1;
  int c, temp;
  
  mMaxScore = -1; //reset
  //look through the last row and find the best score and col
  for(int i=0; i<mCols; i++) {
    if(mScoreTable[r][i]>mMaxScore) {
      mMaxScore=mScoreTable[r][i];
      c=i;
    }    
  }

  if(mMaxScore==-1)
    return;

  //reset gap count
  mSGapCount=0;
  mTGapCount=0;
  mMismatchCount=0;
  mAlignS.clear();
  mAlignT.clear();
  mTAlignEnd=c;

  // Reconstruct the reverse, then flip it. 
  while((r >0) && (c >0)){
	 if (IsSet(DIAG,mPtrTable[r][c])){
		mAlignS.push_back(mS[r-1]);
		mAlignT.push_back(mT[c-1]);
		if(mS[r-1] != mT[c-1]) mMismatchCount++;
		c--;
		r--;
	 }else{
		if (IsSet(TOP,mPtrTable[r][c])){
		  mAlignS.push_back(mS[r-1]);
		  mAlignT.push_back(GAP);
		  r--;
		  mTGapCount++;
		}else{
		  if (IsSet(LEFT,mPtrTable[r][c])){
			 mAlignS.push_back(GAP);
			 mAlignT.push_back(mT[c-1]);
			 c--;
			 mSGapCount++;
		  }
		}
	 }
  }
	
  // Since we fill back to front, when the second strings is not
  // at the beginning,  will be
  // filled with blanks.  This fills in the alignment for 
  // this section. 
  if (r > 0){
	 for(int i = r; i> 0; i--) {
		mAlignS.push_back(mS[i-1]);
		mAlignT.push_back(GAP);
		mTGapCount++;
	 }
  }

  // Reverse, since we constructed a backward version.  
  mAlignS.Reverse();
  mAlignT.Reverse();

  //update the start and end of the alignment which are the start and
  //end of the query and target
  mAlignStart=1;
  mAlignEnd=mS.size();
  mTAlignStart=mTAlignEnd-mAlignT.size()+mTGapCount+1;
}

//------------------------------------------------------
// * RecoverBandInternalAlignment
//------------------------------------------------------
// Reconstruction of alignment from pointers. 
// only accept the alignment that end at the end of the second
// seq, that is, the row=mRow-1
// 
void SequenceAlignment::RecoverBandInternalAlignment()
{
  int r=mBRows-1;
  int c, temp, vrow, vcol;
  
  mMaxScore = -1; //reset
  //look through the last vrow and vcol and find the best score and col
  //if the query stick out, look for vrow and vcol
  if(r==mS.size()) { //the band end touch query end
    for(int i=1; i<mBCols-1; i++) {
      if( (i+r+mTBandStart-2)>mT.size() ) break;
      if(mbandScoreTable[r][i]>mMaxScore) {
	mMaxScore=mbandScoreTable[r][i];
	c=i;
	vrow=mBRows-1;
      }    
    }
  }

  if( (mTBandEnd+mS.size()-1)>mT.size() ) {
    for(int i=mBRows-1; i>0; i--) {
      if( (i-2+mTBandStart+mBCols)<=mT.size() ) break;  //scan for vcol==mT.size done
      for(int j=1; j<mBCols; j++) {
	if( (j+i-2+mTBandStart) == mT.size() ) {
	  if(mbandScoreTable[i][j]>mMaxScore) {
	    mMaxScore=mbandScoreTable[i][j];
	    c=j;
	    vrow=i;
	  }
	}
      }    
    }
  }
  
  if(mMaxScore == -1)
    return;

  //reset gap count
  mSGapCount=0;
  mTGapCount=0;
  mMismatchCount=0;
  mAlignS.clear();
  mAlignT.clear();
  vcol=mTAlignEnd=c+vrow-2+mTBandStart; //the actual col number, start with 1

  //if the query is longer that the target, first push into some gap
  if(vrow<mS.size()) { //march from the end of query to the target end
    for(int i = mS.size(); i>vrow; i--) {
       mAlignS.push_back(mS[i-1]);
       mAlignT.push_back(GAP);
       mTGapCount++;
       mMaxScore += mGapScore;
     }
  }
 
  // Reconstruct the reverse, then flip it.
  r=vrow;
  while((r >0) && (vcol >0)){
    c=vcol-mTBandStart+2-r;
    if (IsSet(DIAG,mbandPtrTable[r][c])){
      mAlignS.push_back(mS[r-1]);
      mAlignT.push_back(mT[vcol-1]);
      if(mS[r-1] != mT[vcol-1]) mMismatchCount++;
      vcol--;
      r--;
    }else{
      if (IsSet(TOP,mbandPtrTable[r][c])){
	mAlignS.push_back(mS[r-1]);
	mAlignT.push_back(GAP);
	r--;
	mTGapCount++;
      }else{
	if (IsSet(LEFT,mbandPtrTable[r][c])){
	  mAlignS.push_back(GAP);
	  mAlignT.push_back(mT[vcol-1]);
	  vcol--;
	  mSGapCount++;
	}
      }
    }
  }
	
  // Since we fill back to front, when the second strings is not
  // at the beginning,  will be
  // filled with blanks.  This fills in the alignment for 
  // this section. 
  if (r > 0){
    for(int i = r; i> 0; i--) {
      mAlignS.push_back(mS[i-1]);
      mAlignT.push_back(GAP);
      mTGapCount++;
    }
  }

  // Reverse, since we constructed a backward version.  
  mAlignS.Reverse();
  mAlignT.Reverse();

  //update the start and end of the alignment which are the start and
  //end of the query and target
  mAlignStart=1;
  mAlignEnd=mS.size();
  mTAlignStart=mTAlignEnd-mAlignT.size()+mTGapCount+1;
}


//------------------------------------------------------
// * RecoverBandEndAlignment
//------------------------------------------------------ 
int SequenceAlignment::RecoverBandEndAlignment()
{
  int r=mBRows-1;
  int c, temp, vrow, vcol;
  
  mMaxScore = -1; //reset
  //look through the last vrow and vcol and find the best score and col
  //if the query stick out, look for vrow and vcol
  if(r==mS.size()) { //the band end touch query end
    for(int i=1; i<mBCols-1; i++) {
      if( (i+r+mTBandStart-2)>mT.size() ) break;
      if(mbandScoreTable[r][i]>mMaxScore) {
	mMaxScore=mbandScoreTable[r][i];
	c=i;
	vrow=mBRows-1;
      }    
    }
  }

  if( (mTBandEnd+mS.size()-1)>mT.size() ) {
    for(int i=mBRows-1; i>0; i--) {
      if( (i-2+mTBandStart+mBCols)<=mT.size() ) break;  //scan for vcol==mT.size done
      for(int j=1; j<mBCols; j++) {
	if( (j+i-2+mTBandStart) == mT.size() ) {
	  if(mbandScoreTable[i][j]>mMaxScore) {
	    mMaxScore=mbandScoreTable[i][j];
	    c=j;
	    vrow=i;
	  }
	}
      }    
    }
  }
  
  if(mMaxScore == -1)
    return 0;

  //reset gap count
  mSGapCount=0;
  mTGapCount=0;
  mMismatchCount=0;
  mAlignS.clear();
  mAlignT.clear();
  vcol=mTAlignEnd=c+vrow-2+mTBandStart; //the actual col number, start with 1

  // Reconstruct the reverse, then flip it.
  r=vrow;
  while((r >0) && (vcol >0)){
    c=vcol-mTBandStart+2-r;
    if (IsSet(DIAG,mbandPtrTable[r][c])){
      mAlignS.push_back(mS[r-1]);
      mAlignT.push_back(mT[vcol-1]);
      if(mS[r-1] != mT[vcol-1]) mMismatchCount++;
      vcol--;
      r--;
    }else{
      if (IsSet(TOP,mbandPtrTable[r][c])){
	mAlignS.push_back(mS[r-1]);
	mAlignT.push_back(GAP);
	r--;
	mTGapCount++;
      }else{
	if (IsSet(LEFT,mbandPtrTable[r][c])){
	  mAlignS.push_back(GAP);
	  mAlignT.push_back(mT[vcol-1]);
	  vcol--;
	  mSGapCount++;
	}
      }
    }
  }
	
  // Reverse, since we constructed a backward version.  
  mAlignS.Reverse();
  mAlignT.Reverse();

  //update the start and end of the alignment which are the start and
  //end of the query and target
  mAlignStart=vrow-mAlignS.size()+mSGapCount+1;
  mAlignEnd=vrow;
  mTAlignStart=mTAlignEnd-mAlignT.size()+mTGapCount+1;

  // Since we fill back to front, when the second strings is not
  // at the beginning,  will be
  // filled with blanks.  This fills in the alignment for 
  // this section. 
  if (vcol > 0 && (vrow == mS.size()) ){
    cout<<"The query is completely internal of the target sequence"<<endl;
    return 0;
  }

  return 1; //yes, query is overlap to the target end
}


//------------------------------------------------------
// * RecoverEndAlignment
//------------------------------------------------------
// Reconstruction of alignment from pointers. 
// only accept the alignment that end at the end of the
// seq, that is, the row=mRow-1 or col=mCol-1
// 
int 
SequenceAlignment::RecoverEndAlignment()
{
  int r = mRows-1, c=mCols-1;
  int init_r, init_c;
  
  mMaxScore = -1; //reset
  //look through the last row and col to find the best score and row and col
  for(int i=0; i<mCols; i++) {
    if(mScoreTable[r][i]>mMaxScore) {
      mMaxScore=mScoreTable[r][i];
      c=i;
    }    
  }
  for(int i=0; i<mRows; i++) {
    if(mScoreTable[i][mCols-1]>mMaxScore) {
      mMaxScore=mScoreTable[i][mCols-1];
      c=mCols-1;
      r=i;
    }    
  }
  //record the init row and col
  init_r=r;
  init_c=c;

  //reset gap count
  mSGapCount=0;
  mTGapCount=0;
  mMismatchCount=0;
  mAlignS.clear();
  mAlignT.clear();
  mTAlignEnd=c;

  // Reconstruct the reverse, then flip it. 
  while((r >0) && (c >0)){
	 if (IsSet(DIAG,mPtrTable[r][c])){
		mAlignS.push_back(mS[r-1]);
		mAlignT.push_back(mT[c-1]);
		if(mS[r-1] != mT[c-1]) mMismatchCount++;
		c--;
		r--;
	 }else{
		if (IsSet(TOP,mPtrTable[r][c])){
		  mAlignS.push_back(mS[r-1]);
		  mAlignT.push_back(GAP);
		  r--;
		  mTGapCount++;
		}else{
		  if (IsSet(LEFT,mPtrTable[r][c])){
			 mAlignS.push_back(GAP);
			 mAlignT.push_back(mT[c-1]);
			 c--;
			 mSGapCount++;
		  }
		}
	 }
  }

  // Reverse, since we constructed a backward version.  
  mAlignS.Reverse();
  mAlignT.Reverse();

  //update the start and end of the alignment which are the start and
  //end of the query and target
  mAlignStart=init_r-mAlignS.size()+mSGapCount+1;
  mAlignEnd=init_r;
  mTAlignStart=mTAlignEnd-mAlignT.size()+mTGapCount+1;

  // Since we fill back to front, when the second strings is not
  // at the beginning,  will be
  // filled with blanks.  This fills in the alignment for 
  // this section. 
  if (c > 0 && (init_r == mRows-1) ){
    cout<<"The query is completely internal of the target sequence"<<endl;
    return 0;
  }

  return 1; //yes, query is overlap to the target end
 
}
//------------------------------------------------------
// * RecoverLocalAlignment
//------------------------------------------------------
// Reconstruction of alignment from pointers. 
//  
// KJD: Assumes that larger string is in mT
// 
void 
SequenceAlignment::RecoverLocalAlignment()
{
	int r = mMaxRow,c = mMaxCol;

	mAlignS.clear();
	mAlignT.clear();
	//reset gap count
	mTGapCount=0;
	mSGapCount=0;
	mMismatchCount=0;

	// The alignment end is with resptect to the longer string, 
	// here assumed to be the one that corresponds to columns. 
	mAlignEnd = mMaxRow;
 
	while((r > 0) && (c > 0) && (mScoreTable[r][c]!=0)){
	  if (IsSet(DIAG,mPtrTable[r][c])){
		 mAlignS.push_back(mS[r-1]);
		 mAlignT.push_back(mT[c-1]);
		 if(mS[r-1] != mT[c-1]) mMismatchCount++;
		 c--;
		 r--;
	  }else{
		 if (IsSet(TOP,mPtrTable[r][c])){
			mAlignS.push_back(mS[r-1]);
			mAlignT.push_back(GAP);
			mTGapCount++;
			r--;
		 }else{
			if (IsSet(LEFT,mPtrTable[r][c])){
			  mAlignS.push_back(GAP);
			  mAlignT.push_back(mT[c-1]);
			  mSGapCount++;
			  c--;
			}
		 }
	  }		
	}

	//if (mAlignS.size() < mS.size()){
	//  cout<<"Alignment Size Mismatch :"<<r<<","<<c<<endl;
	//  mScoreTable.Print(r-1,r+10,c-1,c+10);
	//}else{
	//  cout<<"Alignment Size OK :"<<r<<","<<c<<endl;
	//  mScoreTable.Print(r-1,r+10,c-1,c+10);
	//}

	// The reverse has been constructed, via push_back, so we 
	// flip it to get the original sense. 
	mAlignS.Reverse();
	mAlignT.Reverse();
	
	// Compute the start of the alignment, which is the end minus
	// the size of the shorter string. 
	mAlignStart = mAlignEnd - mAlignS.size()+mSGapCount+1;
	mTAlignEnd = mMaxCol;
	mTAlignStart = mTAlignEnd - mAlignT.size()+mTGapCount+1;
}


//------------------------------------------------------
// * PrintSpacePaddedNum
//------------------------------------------------------
//
void
SequenceAlignment::PrintSpacePaddedNum(int x)
{
	// Add a space to all positive numbers to match negative nums. 
	if (x>=0){
		cout<<" ";
	}

	if (abs(x) < 10){
		cout<<"   ";
	}else{
		if (abs(x) < 100){
			cout<<"  ";
		}else{
			if (abs(x) < 1000){
				cout<<" ";
			}
		}
	}
	cout<<x;
}

//------------------------------------------------------
// * PrintScoreTable
//------------------------------------------------------
//
//
void
SequenceAlignment::PrintScoreTable()
{
	int i,j;
	
	// print heading. 
	cout<<"            ";
	for(j = 0;j<mT.size();j++){
		cout<<mT[j]<<"     ";
	}
	cout<<"\n";
	
	for(i = 0;i<mRows;i++){
		if (i > 0){
			cout<<mS[i-1]<<":";
		}else{
			cout<<"  ";
		}
		for(j = 0;j<mCols;j++){
			PrintSpacePaddedNum(mScoreTable[i][j]);
			cout<<" ";
		}
		cout<<"\n";
	}
}

//------------------------------------------------------
// * PrintbandScoreTable
//------------------------------------------------------
//
//
void
SequenceAlignment::PrintBandScoreTable()
{
	int i,j;
	
	// print heading. 
	cout<<"            ";
	for(j = 0;j<mT.size();j++){
		cout<<mT[j]<<"     ";
	}
	cout<<"\n";
	
	for(i = 0;i<mBRows;i++){
		if (i > 0){
			cout<<mS[i-1]<<":";
		}else{
			cout<<"  ";
		}
		for(j = 0;j<mBCols;j++){
			PrintSpacePaddedNum(mbandScoreTable[i][j]);
			cout<<" ";
		}
		cout<<"\n";
	}
}

//------------------------------------------------------
// * PrintPtrTable
//------------------------------------------------------
//
void
SequenceAlignment::PrintPtrTable()
{
	int i,j;
	
	// print heading. 
	cout<<"            ";
	for(j = 0;j<mT.size();j++){
		cout<<mT[j]<<"      ";
	}
	cout<<"\n";

	for(i = 0;i<mRows;i++){
		if (i > 0){
			cout<<mS[i-1]<<":";
		}else{
			cout<<"  ";
		}
		cout<<"   ";
		for(j = 0;j<mCols;j++){
			if (IsSet(LEFT,mPtrTable[i][j])) cout<<"L";
			else cout<<".";
				
			if (IsSet(DIAG,mPtrTable[i][j])) cout << "D";
			else cout<<".";
			
			if (IsSet(TOP,mPtrTable[i][j])) cout<<"T   ";
			else cout<<".   ";
			cout<<" ";
		}
		cout<<"\n";
	}
}
	

//------------------------------------------------------
// * Clear
//------------------------------------------------------
// 
void
SequenceAlignment::Clear(int signal,int &c)
{
  c = c & (~signal);
}

//------------------------------------------------------
// * IsSet
//------------------------------------------------------
// 
inline bool
SequenceAlignment::IsSet(int signal,int c)
{
  return(c & signal);
}

//------------------------------------------------------
// * IsSet
//------------------------------------------------------
// 
bool
SequenceAlignment::IsSet(int signal,char c)
{
  return(c & signal);
}


//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
inline void
SequenceAlignment::Set(int signal,char &c)
{
  c = c | signal;
}

//------------------------------------------------------
// * Set
//------------------------------------------------------
// 
inline void
SequenceAlignment::Set(int signal,int &c)
{
  c = c | signal;
}

//------------------------------------------------------
// * Print
//------------------------------------------------------
//
// 
void
SequenceAlignment::Print(int start,int end)
{
	int i;
	for(i = start;i<end;i++){
		cout<<mAlignS[i];
	}
	cout<<"\n";
	for(i = start;i<end;i++){
		cout<<mAlignT[i];
	}
	cout<<"\n";
}

//------------------------------------------------------
// * Print
//------------------------------------------------------
// 
void
SequenceAlignment::Print()
{
	int i;

	for(i = 0;i<mAlignT.size();i++){
		cout<<mAlignT[i];
	}
	cout<<"\n";
	for(i = 0;i<mAlignS.size();i++){
		cout<<mAlignS[i];
	}
	cout<<"\n";
}

//--------------------------------------------------
// * PrintAlign the alignment of two sequence, add on by rui
//----------------------------------------------
//
void SequenceAlignment::PrintAlign(ostream& output) {
  int LongQuery;
  int identity;

  string querystr, targetstr;
  int temp, current; //var for find the size of index for printing
  int templen=0;

  int seqperline, total, i;

  output<<">"<<mT.mName<<mT.mDescription<<"\n";
  output<<"\tLength = "<<mT.size()<<"\n\n";
  
  identity=mAlignS.size()-mMismatchCount-mSGapCount;
  output<<" Score="<<mMaxScore<<", Identity="<<identity<<"/"<<mAlignS.size()
        <<", gap="<<mSGapCount<<"\n\n"; //get bits score later

  //alignment start
  //get the bigger index
  temp = Max(mAlignStart, mAlignEnd, mTAlignStart, mTAlignEnd);

  //find out the length of temp
  for(current=1; current<temp; current=current*10) {
    templen++;
  }
  
  seqperline = 71 - templen*2; //number of nt to put per line
  total=ceil( (double) mAlignS.size()/seqperline);
  output.setf( ios::left, ios::adjustfield);
  mAlignS.CopyToString(querystr);
  mAlignT.CopyToString(targetstr);
  int CurrMStart=mAlignStart;
  int CurrTStart=mTAlignStart;
  int MGap, TGap;
  for(current=0; current<total; current++) {
    int trueseqperline=seqperline;
    if( (current*seqperline+seqperline)> mAlignS.size() ) trueseqperline=mAlignS.size()-current*seqperline;
     //a loop to get the count of gaps
    MGap=TGap=0;
    for(i=0; i<trueseqperline; i++) {
      if(mAlignS[i+current*seqperline] == GAP)
	MGap++;
      if(mAlignT[i+current*seqperline] == GAP) 
	TGap++;
    }
    output<<"Query: "<< setw(templen)<<CurrMStart<<" "<<querystr.substr(current*seqperline, trueseqperline);
    output<<" "<<(CurrMStart+trueseqperline-1-MGap)<<"\n";
    CurrMStart=CurrMStart+trueseqperline-MGap;
    
    output<<"        "<<setw(templen)<<" ";
    //a loop to check if there is a mismatch
    for(i=0; i<seqperline; i++) {
      if( (i+current*seqperline) >= mAlignS.size()) break;  //break if reach the end

      if(mAlignS[i+current*seqperline] == mAlignT[i+current*seqperline]) {
	output<<"|";
      } else {
	output<<" ";
      }
    }
    
    output<<"\nTargt: "<<setw(templen)<<CurrTStart<<" "<<targetstr.substr(current*seqperline, trueseqperline);
    output<<" "<<(CurrTStart+trueseqperline-1-TGap)<<"\n\n";
    CurrTStart=CurrTStart+trueseqperline-TGap;
  } //end of output for current
  output<<"\n\n";

}

//------------------------------------------------------
// * Max
//------------------------------------------------------
// 
inline int
SequenceAlignment::Max(int a,int b,int c)
{
  if (a>b){
	 if (a>c){
		return(a);
	 }
  }else{ // a < b
	 if (b>c){
		return(b);
	 }else{ // a<b b<c
		return(c);
	 }
  }
  return(c);
}


//------------------------------------------------------
// * Max
//------------------------------------------------------
// 
inline int
SequenceAlignment::Max(int d,int b,int c,int a)
{
  if (a>b){
	 if (a>c){
		if(a>d){
		  return(a);
		}else{
		  return(d);
		}
	 }else{ // a>b a<c 
		if (c>d){
		  return(c);
		}else{
		  return(d);
		}
	 }
  }else{ // a<b
	 if (b>c){
		if (b>d){
		  return(b);
		}else{
		  return(d);
		}
	 }else{ // a<b b<c
		if (c>d){
		  return(c);
		}else{
		  return(d);
		}
	 }
  }
}

//Score fun
int SequenceAlignment::Score(char A, char B)
{
    if (A == B){
      return mMatchScore;
    } else if( IgnorN() ) {
      return mMismatchScore;
    } else if( A == 'N' || A == 'n' || B == 'N' || B == 'n') { 
      return mMatchScore;
    }
    return mMismatchScore;
}














